from pylab import *
import scipy.optimize as sp
from pitubalib import *

# Parameters to optimize:
# Bo
# Rs
# Pb
# VOIP


#---------------------------------------------------------------------------
# PVT FROM SERGI RESERVOIR SAMPLE
#---------------------------------------------------------------------------

## PVT
## PVT do Bloco 3
fluido=classPVT()
fluido.pe=array([200.,170.,140.,110.,80.,60.,48.,35.,25.,15.,1.033])
fluido.Boe=0. # PARAM
fluido.Rse=0. # array([66.9,57.8,48.6,39.4,30.2,24.1,20.4,16.4,13.1,9.4,0.0000])
fluido.Bge= array([0.0047,0.055,0.0067,0.0085,0.0117,0.0156,0.0195,0.0310,0.0435,0.0720,1.0000])
fluido.Bwe=1.0*ones(len(fluido.pe))
fluido.viscoe=array([5.81000,7.23000,8.05000,9.05000,10.23000,11.58000,12.17000,12.80000,14.12000,14.83000,15.57000,16.33000,17.12000,17.96000,18])
fluido.viscwe=1.0*ones(len(fluido.pe))
fluido.co=9.4e-5
fluido.Pb=0. #48.
fluido.T=273.+57.


#---------------------------------------------------------------------------
# COMPUTING EQUATION COEFFICIENTS
# Np*(Bo2+(Rp-Rso2)*Bg2)+Wp*Bw2=N*(Eo+m*Eg+Efw)+(Wi+We)*Bw2+Gi*Bg2
#---------------------------------------------------------------------------
Eq=classEBM()
prod=genfromtxt('prod_AG_BL4.txt',delimiter='\t',skiprows=1)
Eq.Np=prod[:,0]*1e3
Eq.Gp=0.1*prod[:,1]*1e3
Eq.Wp=prod[:,2]*1e3
Eq.Wi=prod[:,3]*1e3
Eq.dt=prod[:,4]
Eq.We=Eq.Wp-Eq.Wi

for i in range(1,len(Eq.We)):
	Eq.We[i]=max(Eq.We[i],Eq.We[i-1])

Eq.N=0.75e6

AG=classRes()
# DEFINING DATUM=-790M (AVERAGE RESERVOIR DEPTH)
# Pasta do poco CER-4: pe=91.3 kgf/cm2 @ MD=940m, COTA=-835.4,
AG.p0= 100. #91.3+0.08*(-835.4+800.)
AG.Pb0=fluido.Pb

AG.phi=0.21
AG.Sw0=0.13
AG.Swi=0.13
AG.Sor=0.25



# Petrophysics
AG.m=0.
AG.cr=2.5e-5 #3.5e-6 * 14.7/1.033#5e-5
petroPhys=classPetrophys()
Swe=array([0.13,0.23,0.34,0.44,0.54,0.65,0.75])
petroPhys.Swe=Swe
petroPhys.kroe=0.7*((1.-AG.Sor-Swe)/(1.-AG.Sw0-AG.Sor))**3
petroPhys.krwe=0.075*((Swe-AG.Sw0)/(1.-AG.Sw0-AG.Sor))**6

# Aquifer
aquifero=classAquifero()
aquifero.model='user_Np'
aquifero.Npe=Eq.Np
aquifero.Wee=Eq.We
aquifero.WW=15*Eq.N # 400e6*0.27 # 35e6
aquifero.p0=AG.p0
aquifero.cr=AG.cr
aquifero.cw=5.7e-5
aquifero.pvt=fluido
aquifero.k=200.
aquifero.L=1000.


# TESTES DE PRESSAO

# DATUM=-800m


AG.teste.Np=array([0.814,4.35,21.64,36.4095,38.9854,78.1342,95.2744,108.4334,128.865,143.653,164.3191,171.2098,194.5308])*1e3
AG.teste.p=array([91.5,70.41,54,48.89,48.685,32.6,25.69,20.2,13.085,10.1,17.99,19.485,22.855])
##AG.teste.p=array([72.69,53.12,47.73,47.4,28.71,24.52,21.29,11.55,11.08,16.61,16.28,20.5])*1.
##AG.teste.Np=array([3.87430,24.67000,36.40000,42.44000,65.32000,95.85000,108.40000,129.38000,144.10000,164.35000,171.40000,194.5])*1e3

#AG.teste.Gp=interp(AG.teste.Np,Eq.Np,Eq.Gp)
AG.teste.Wp=interp(AG.teste.Np,Eq.Np,Eq.Wp)
AG.teste.Wi=interp(AG.teste.Np,Eq.Np,Eq.Wi)

Eq.res=AG
Eq.aquifero=aquifero
Eq.pvt=fluido
Eq.petroPhys=petroPhys


#-----------------------------------------------------------------------
# Optimize parameters
#-----------------------------------------------------------------------
# COMPUTE GAS PRODUCTION FROM EXP PRESSURE

def trial(x,Eq):
	print '-------------------------'
	print x

	Bo=x[0]
	Rs=x[1]
	Pb=x[2]
	N=x[3]
	co=x[4]
	n=len(Eq.pvt.pe)
	Eq.pvt.Boe=linspace(Bo,1.,n)
	Eq.pvt.Rse=Rs/(Pb-Eq.pvt.pe[-1])*(Eq.pvt.pe-Eq.pvt.pe[-1])
	Eq.pvt.Pb=Pb
	Eq.pvt.co=co
	Eq.res.Pb0=Pb
	Eq.res.N=N
	Eq.res.Vp0=Eq.N*Eq.pvt.Bo(Eq.res.p0,Eq.res.Pb0)/(1.-Eq.res.Sw0-Eq.res.Sg0)
	Gp=Eq.calcGp()
	Eq.Gp=interp(Eq.Np,Gp[:,0],Gp[:,1])
	savetxt('Gp.dat',Gp)
	out=Eq.runHist(1.,Eq.res.p0)

	residual=sum((Eq.res.teste.p-interp(Eq.res.teste.Np,Eq.Np,array(out.p)))**2)

	print 'Residual =',residual
	print len(out.p),len(Eq.Np)
	print '-------------------------'
	return residual

def approx_fprime(xk,Eq):
	epsilon=0.02*xk
	f0 = trial(xk,Eq)
	grad = zeros(len(xk))
	ei = zeros(len(xk))
	for k in range(len(xk)):
		ei[k]=epsilon[k]
		fu=trial(xk+ei,Eq)
		fd=trial(xk-ei,Eq)
		grad[k] = (fu - fd)/ei[k]/2.
		ei[k] = 0.0
	return grad


x0=array([1.2,80.,60.,0.75e6,1e-5])

##xMin= sp.brute(trial, ((1.05,1.25),(20.,60.),(40.,70.),(0.7e6,1.2e6)), args=[Eq], Ns=5, full_output=0)
##xMin= sp.fmin_cg(trial, x0, fprime=approx_fprime, args=[Eq],epsilon=0., gtol=0.1, norm=inf, maxiter=None, full_output=0, disp=1, retall=0, callback=None)
xMin = sp.fmin_l_bfgs_b(trial, x0, fprime=approx_fprime, args=[Eq], approx_grad=0, bounds=((1.05,1.2),(10.,25.),(20.,65.),(0.6e6,1.2e6),(5e-5,1.5e-4)), m=10, factr=10000000.0, pgtol=1e-05, epsilon=1e-08, iprint=-1, maxfun=15000, disp=None)
xMin=xMin[0]

#xMin=array([1.25000000e+00,2.50000000e+01,5.68683375e+01,7.50000000e+05,1.75205742e-04])
print 'xMin = ',xMin


Bo=xMin[0]
Rs=xMin[1]
Pb=xMin[2]
N=xMin[3]
co=xMin[4]

n=len(Eq.pvt.pe)
Eq.pvt.Boe=linspace(Bo,1.,n)
Eq.pvt.Rse=Rs/(Pb-Eq.pvt.pe[-1])*(Eq.pvt.pe-Eq.pvt.pe[-1])
Eq.pvt.Pb=Pb
Eq.pvt.co=co
Eq.res.Pb0=Pb
Eq.res.N=N

print 'Rs = ',Eq.pvt.Rse[0]
Eq.res.N=N
Eq.pvt.co=co
Eq.res.Vp0=Eq.N*Eq.pvt.Bo(Eq.res.p0,Eq.res.Pb0)/(1.-Eq.res.Sw0-Eq.res.Sg0)
Gp=Eq.calcGp()
Eq.Gp=interp(Eq.Np,Gp[:,0],Gp[:,1])
savetxt('Gp.dat',Gp)
out=Eq.runHist(1.,Eq.res.p0)

figure()
plot(Eq.Np*1e-6,out.p,'bo-')
plot(Eq.Np*1e-6,out.Pb,'--k')
errorbar(AG.teste.Np*1e-6,AG.teste.p,fmt='ro',yerr=60./14.23)
title('Material Balance :: CER Oil Field, Block: 4, Zone: AG.')
legend(('Calc Press','Bubble Press','Measured'))
xlabel('Cumulative oil production (1e6 std m3)')
ylabel('Reservoir Pressure (kgf/cm2)')
ylim([0,100])
grid(1)

figure()
plot(Eq.Np*1e-6,gradient(Eq.Gp)/gradient(Eq.Np),'rs-')
plot(Eq.Np*1e-6,zeros(len(Eq.Np)),'k:')
xlabel('Np (1e6 std m3)')
ylabel('GOR')
grid(1)

figure()
plot(Eq.Np*1e-6,Eq.Gp*1e-6,'bo-')
title('Material Balance :: CER Oil Field, Block: 4, Zone: Agua Grande.')
xlabel('Cumulative oil production (1e6 std m3)')
ylabel('Cumulative Gas production (1e6 std m3)')

figure()
plot(Eq.Np*1e-6,array(Eq.Wi)*1e-6,'bx-')
plot(Eq.Np*1e-6,array(out.We)*1e-6,'gx-')
title('Material Balance :: CER Oil Field, Block: 4, Zone: Agua Grande.')
legend(('Wi','We'),'upper left')
xlabel('Cumulative oil production (1e6 std m3)')
ylabel('Water Volume (1e6 std m3)')

figure()
plot(Eq.Np*1e-6,out.Sw,'bo-')
plot(Eq.Np*1e-6,out.So,'ko-')
plot(Eq.Np*1e-6,out.Sg,'ro-')
plot(Eq.Np*1e-6,array(out.Sg)+array(out.So)+array(out.Sw),'k-')
xlabel('Cumulative oil production (1e6 std m3)')
legend(('Sw','So','Sg','Sw+So+Sg'),'center left')
grid(1)
title('Saturations')

#figure()
#semilogy(range(len(Eq.Np)),abs(array(out.residual)),'k-') # Normalization due to sweep efficiency
#xlabel('step')
#ylabel('Residual')
#grid(1)

figure()
plot(1.-array(out.So),out.fw,'ro') # Normalization due to sweep efficiency
xlabel('1-So')
ylabel('BSW')
#xlim([0.,1.])
grid(1)

# Computing reservoir IP
#PI=gradient(Eq.Wp+Eq.Np)/array(out.p)*1e-3
#figure()
#plot(Eq.Np*1e-6,PI,'k-')
#xlabel('Cumulative oil production (1e6 std m3)')
#legend(('Reservoir PI'),'center left')
#grid(1)
#title('Reservoir PI')


#-----------------------------------------------------------------------
show()
##print 'Gp max = ',(Eq.N*(fluido.Rs(AG.Pb0,AG.Pb0)-fluido.Rs(20.,AG.Pb0))+Eq.Np[-1]*fluido.Rs(20.,AG.Pb0))/1e6
##print 'Porous Volume = ',AG.Vp0/1e6
##print 'Savg = ',AG.Vp0/18.
##print 'hgas = ',AG.Vp0*out.Sg[-1]/18.
